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EXECUTIVE  SUMMARY 


This  work  investigated  the  detection  of  pulses  and  extraction  of  modulation 
parameters  from  different  types  of  time-limited  chirp  signals,  as  may  be  found  in  pulse- 
compression  radar  signals.  The  work  is  split  into  two  parts.  The  first  part  examines  the 
pulse  detection  problem,  i.e.,  the  detection  of  the  pulse  start/stop  times.  Such  information 
can  be  used  to  determine  the  pulse  width  and  repetition  rate  of  the  radar  systems  under 
investigation  in  an  automated  fashion.  We  compare  the  robustness  of  three  TCF-based 
schemes  and  an  envelope  detection  algorithm  in  noisy  environments.  Results  show  that 
none  of  the  pulse  detection  schemes  considered  in  this  work  to  be  clearly  better  than  the 
others  for  the  SNR  range  considered,  and  the  specific  selection  to  be  a  function  of  the 
desirable  characteristic,  either  Pfa  or  Pd,  to  be  optimized. 

The  second  part  of  the  work  focuses  on  the  extraction  of  modulation  parameters 
from  two  specific  modulation  types:  linear  and  hyperbolic  chirp  modulation.  For  the 
second  part  it  is  assumed  that:  1)  individual  pulses  have  been  detected  and  isolated  prior 
to  the  processing,  and  2)  the  modulation  type  is  known,  as  we  do  not  discuss  issues 
related  to  the  specific  identification  of  the  modulation  scheme  which  are  left  for  further 
research.  The  main  idea  behind  the  work  reported  in  this  section  is  the  fact  that  we  use  an 
image  processing  approach.  First,  time-frequency/scale  images  of  the  extracted  pulses  are 
generated,  next  additional  processing  applied  to  the  images  allow  the  identification  of  the 
modulation  parameters.  Thus,  different  types  of  time-frequency  and  wavelet-based 
transformations  are  considered  and  applied  to  linear  chirp  signals  in  additive  white 
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Gaussian  noise.  Results  show  that  time-frequency  transformations  lead  to  better  focused 
images  when  dealing  with  noisy  signals,  and  thereby  to  better  estimation  of  the 
modulation  parameters  than  the  wavelet-based  decompositions  do.  Simulations 
investigate  the  robustness  of  the  modulation  estimation  schemes  in  noisy  environments 
both  for  linear  and  hyperbolic  chirps. 


I.  INTRODUCTION 


A.  BACKGROUND 

Numerous  schemes  have  been  proposed  over  the  years  to  retrieve  modulation 
parameters  of  chirp  signals  for  various  types  of  applications  ranging  from 
communications  to  radar  and  sonar.  The  capability  to  extract  the  modulation  parameters 
in  an  automated  fashion  can  be  very  useful  as  such  procedures  can  be  integrated  in 
monitoring  schemes.  The  report  focuses  on  signals  with  constant  amplitudes  and  various 
types  of  modulation  types,  such  as  constant,  linear  and  hyperbolic  intra-pulse  frequency 
modulation  over  a  given  pulse.  These  modulation  types  can  be  found  in  pulse- 
compression  radar  signals,  pulse  or  CW  radar  signals  received  from  targets,  or 
modulation  signals  found  in  communication  applications.  Radar  intra-pulse  modulation 
parameters  are  usually  extracted  using  hardware  schemes,  as  they  are  well  suited  to 
extract  the  instantaneous  frequency  at  high  SNR  levels.  This  study  considers  the  problem 
from  a  different  angle  and  investigates  the  application  of  temporal  correlation  functions 
(TCF),  time-frequency  (TF)  and  time-scale  (wavelet)  transformations,  and  basic  image 
processing  techniques  to  detect  pulse  location  and  extract  modulation  information. 

B.  STUDY  ORGANIZATION 

This  work  was  split  into  two  parts.  The  first  part  examined  the  pulse  detection 
problem,  i.e.,  the  detection  of  the  pulse  start/stop  times,  and  compared  several  TCF-based 
detection  schemes.  This  information  can  be  used  to  determine  the  pulse  width  and 
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repetition  rate  of  the  radar  systems  under  investigation.  The  second  part  of  the  work 
focused  on  the  extraction  of  modulation  parameters  from  two  specific  modulation  types: 
linear  and  hyperbolic  chirp  modulation.  For  the  second  part  it  is  assumed  that:  1) 
individual  pulses  have  been  detected  and  isolated  prior  to  the  processing  discussed  next, 
and  2)  the  modulation  type  is  known,  as  we  do  not  discuss  issues  related  to  the  specific 
identification  of  the  modulation  scheme  which  is  left  for  further  research. 

Chapter  II  presents  a  brief  review  of  RADAR  concepts.  Chapter  III  presents  the 
concept  of  temporal  correlation  in  the  context  of  pulse  detection.  Chapter  IV  describes 
the  various  TCF-based  signal  analysis  techniques  considered  and  the  detection 
performances  obtained  for  pulses  distorted  by  white  Gaussian  noise.  Chapter  V  briefly 
describes  the  time-frequency  transformations  considered  in  this  work  to  extract  intra¬ 
pulse  modulation  parameters.  The  application  of  TF  and  Radon  transforms  to  extract  the 
modulation  parameters  from  linear  chirps  is  discussed  in  Chapter  VI.  Next,  we 
considered  the  extraction  of  modulation  parameters  for  hyperbolic  chirps.  We  restricted 
our  study  to  the  top  three  energy  TF  transformations  leading  to  the  best  image  quality  for 
the  linear  chirp,  and  derived  the  estimation  procedure  presented  in  Chapter  VII.  Finally, 
conclusions  and  recommendations  for  further  research  are  presented  in  Chapter  VDI. 
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II.  MODERN  RADAR 


Radar  signal  identification  and  parameter  identification  would  be  simple  if  radars 
operated  at  a  fixed  frequency  with  constant  parameters.  Modem  techniques  have 
improved  the  radar’s  unambiguous  detection  range,  range  resolution,  and  the  radar’s 
susceptibility  to  jamming  while  lowering  the  power  requirements.  These  modifications 
tend  to  make  the  radar  recognition  easier,  but  they  tend  to  make  the  detection  and 
extraction  of  the  parameters  more  difficult. 

The  modifications  tend  to  fall  into  three  general  categories:  RF  agility,  inter-pulse 
and  intra-pulse  modulation.  Some  examples  of  the  techniques  and  what  improvements 
they  can  provide  are  discussed  next. 

A.  RF  AGILITY 

Radiating  frequency  (RF)  agility  is  the  ability  to  change  the  operating  frequency 
on  a  pulse  to  pulse  basis.  There  are  many  advantages  to  using  an  RF  agile  radar,  although 
the  Electronic  Counter  Counter  Measure  (ECCM)  capability  is  generally  the  most 
emphasized  improvement. 

Pulse  frequency  agile  radars,  particularly  when  the  frequencies  are  randomized 
over  the  operating  range  of  the  radar,  are  difficult  to  jam  [42].  RF  agility  is  also  used  to 
reduce  the  angular  error  in  the  radar  caused  by  glint.  Glint  is  produced  when  echoes  are 
treated  as  if  they  are  coming  from  one  scattering  center,  instead  of  an  echo  vector  sum 
returned  from  the  different  scattering  centers  of  an  extended  target.  By  using  multiple 
frequencies,  selecting  the  frequency  with  the  strongest  return  tends  to  eliminate  the  large 
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angle  errors.  The  minimum  frequency  difference  between  pulses  required  for  glint 
reduction  is  given  by 


Pulse  to  pulse  frequency  agility  is  a  successful  technique  to  counter  ECCM.  Radars 
can  also  employ  a  dwell  and  switch  RF  hopping  scheme.  In  this  case,  a  given  frequency 
for  a  number  of  pulses  or  scans  is  used  followed  by  a  jump  to  another  frequency.  The 
frequencies  are  usually  discretely  spaced  in  a  given  band.  The  frequency  hop  can  occur 
after  a  set  number  of  pulses,  after  a  set  scan  time,  or  with  a  PRF  jump. 

B.  INTERPULSE  MODULATION 

Inter-pulse  modulation  refers  to  the  modulation  of  the  PRF  or  PRI  of  a  radar  pulse 
train.  This  combines  the  advantages  of  a  low  with  that  of  a  high  PRF  radar.  These 
techniques  also  minimize  range  eclipsing  which  occurs  when  the  target’s  echo  arrives 
during  the  next  pulse  receive  cycle.  PRF  modulation  does  not  require  complex  circuits,  so 
it  is  often  used  to  improve  the  radar  performance.  Four  common  techniques  are  described 
below  [41]. 

1.  PRF  Switching 

A  PRF  switching  radar  switches  between  a  few  discrete  PRF’s.  Typically  one 
PRF  will  be  used  for  a  set  time  interval  such  as  one  scan  interval  or  a  given  number  of 
pulses,  which  is  then  switched  to  another  PRF  for  the  next  period.  Some  RF  agile  radars 
can  also  hop  in  frequency  as  the  PRF  switch  occurs. 
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Careful  selection  of  the  PRF’s  can  extend  the  blind  speeds  beyond  the  range  of 
expected  values.  Only  two  or  three  discrete  PRF  values  are  necessary  to  resolve  most 
ambiguous  range  problems  [41]. 

2.  Staggered  PRF 

The  PRF  switches  on  a  pulse  to  pulse  basis.  This  technique  is  mainly  used  to 
eliminate  blind  speeds  and  blind  ranges  [41].  Second  time  around  clutter  echoes  do  not 
cancel  since  the  clutter  does  not  appear  at  the  same  range  from  pulse  to  pulse.  It  takes 
several  looks  to  determine  whether  or  not  a  target  is  in  motion  [42], 

3.  Sliding  PRF 

The  PRF  is  continuously  varied,  typically  starting  at  a  low  and  increasing  to  a  high 
PRF  value.  Then  the  PRF  jumps  back  to  the  low  value  to  begin  sliding  again.  Though  the 
low-to-high  variations  tend  to  be  repeated,  the  PRF  values  are  not  always  the  same  [41], 

By  continuously  increasing  the  PRF  the  target  echo  can  not  catch  up  with  the 
range,  and  range  eclipsing  will  not  occur.  Because  of  this.  Sliding  PRF  is  often  used  in 
target  tracking  radar  [41]. 

4.  Jittered  PRF 

The  start  time  of  each  successive  pulse  is  varied  (jittered)  relative  to  where  it 
should  have  started  in  a  regular  pulse  train.  Jitter  is  normally  classified  by  the  percent 
jitter,  which  is  the  maximum  jitter  time  divided  by  the  normal  pulse  interval,  and  by  the 
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manner  in  which  the  jitter  time  varies.  The  jitter  may  vary  randomly  from  pulse  to  pulse 
but  the  actual  variations  may  be  from  a  pre-determined  discrete  set.  The  variation  may 
also  follow  a  sinusoidal,  a  triangular,  or  some  other  prescribed  pattern  [41]. 

Jittering  the  PRF  improves  the  radar’s  performance  both  with  respect  to  blind 
speed  and  the  unambiguous  range.  Some  commonly  used  jitter  time  distributions  can 
improve  both  parameters  by  a  factor  of  two  or  more  [41]. 

C.  INTRA-PULSE  MODULATION 

Intra-pulse  modulation  in  conjunction  with  pulse  compression  delivers  the  power 
of  a  long  pulse  (a  better  maximum  range)  while  keeping  the  range  resolution  of  a  short 
pulse.  There  are  two  methods  to  obtain  intra-pulse  modulation,  frequency  modulation  of 
the  pulse  (FMOP)  and  phase  modulation  of  the  pulse  (PMOP)  [41]. 

Pulse  compressed  radars  have  a  greater  effective  range  and  range  resolution.  They 
also  are  less  susceptible  to  jamming  as  the  jammer’s  power  is  spread  over  a  larger 
bandwidth.  Pulse  compression  radars  are  more  tolerant  of  other  radars  operating  in  the 
same  frequency  spectrum  (mutual  interference).  They  can  operate  in  a  given  spectral 
band  with  another  pulse-compressed  radar  if  each  of  them  has  its  own  characteristic 
modulation  [42]. 

Some  disadvantages  of  intra-pulse  modulation  are  the  increase  in  cost  and 
complexity  of  the  system.  Also,  the  time  sidelobes  arising  from  compression  can  mask  or 
be  mistaken  for  a  target  [41]. 
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Relative  to  short  pulse  radars,  pulse  compressed  radars  have  the  disadvantage  of  a 
longer  minimum  range.  If  the  sidelobes  are  limited  small  signals  may  be  lost.  Finally, 
pulse  compressed  radars  do  not  have  the  short  pulse  radars  natural  immunity  to  repeat 
jammers  and  range  gate  stealing  jammers  [42]. 

1.  Frequency  Modulation  on  the  Pulse  (FMOP) 

A  common  type  of  FMOP  uses  linear  frequency  modulation.  This  modulation, 
also  known  as  a  “chirp”,  linearly  increases  (or  decreases)  the  frequency  of  the  pulse  over 
a  small  frequency  range  while  holding  the  amplitude  and  pulse  width  constant.  The 
receiver  is  designed  so  that  all  the  frequencies  arrive  at  the  detector  at  the  same  time. 
This  compresses  the  pulse  to  a  bandwidth  equal  to  the  frequency  shift  of  the  pulse  [41]. 

Another  prominent  FMOP  technique  is  FM  Stepping.  The  pulse  is  broken  up  into 
sub  pulses  that  are  held  at  a  constant  frequency,  but  stepped  up  or  down  by  a  discrete 
frequency  value.  The  bandwidth  of  the  compressed  signal  is  the  number  of  steps  times 
the  change  in  frequency  [41]. 

2.  Phase  Modulation  on  the  Pulse  (PMOP) 

The  pulse  is  broken  into  N  equal  length  sub  pulses.  Each  of  these  sub  pulses  is 
either  in  phase  or  180  degrees  out  of  phase.  The  phase  of  the  sub  pulses  is  determined  by 
a  preset  code,  usually  a  Barker  code  or  a  pseudo  random  noise  sequence.  A  matched  filter 
cross  correlates  the  received  signal  with  a  template.  When  the  received  signal  and  the 
template  match  up,  the  receiver  returns  a  large  peak  [41], 
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3.  FMOP  and  PMOP  Comparison 

The  PMOP  radar  has  a  poor  resolution  in  a  dense  target  environment.  PMOP 
radars  benefit  from  the  different  phase  coded  sequences  that  can  be  assigned  to  different 
radars,  so  that  a  number  of  PMOP  radars  can  operate  in  the  same  frequency  band.  PMOP 
radars  are  less  susceptible  to  repeat  jammers  than  FMOP  radars  are,  since  the  code  can  be 
changed  to  counter  a  repeater  mimicking  it.  FMOP  radars  are  likely  to  be  used  when  a 
wide  bandwidth  or  very  narrow  compressed  pulse  is  desired.  A  PMOP  radar  is  likely  to 
be  used  when  jamming  is  a  problem  or  when  there  will  be  multiple  radars  using  the  same 
portion  of  the  frequency  band  [42]. 

D.  COMBINING  TECHNIQUES 

While  each  of  the  three  techniques  have  been  discussed  separately  it  is  not 
unusual  for  a  radar  system  to  use  a  combination  of  these  techniques.  It  is  also  important 
to  note  that  a  radar  system  may  have  multiple  operational  modes  which  can  use  any  one 
or  a  combination  of  the  techniques. 

A  typical  example  is  a  missile  system  radar.  It  may  use  a  constant  PRF  while  in  a 
search  mode,  switch  to  an  RF  dwell  and  switch  mode  to  acquire  the  target,  and  then 
switch  to  an  RF  agile,  PRF  staggered,  chirped  pulse  modulation  to  track  a  target  before 
launching  a  missile. 


III.  TEMPORAL  CORRELATION  FUNCTION 

The  temporal  correlation  function  allows  for  non-stationary  situations.  It  is 
indexed  in  time  and  in  delay. 

A.  CORRELATION  FUNCTION  DEFINITION 

Depending  on  the  underlying  process,  various  definitions  are  given  to  the 
auto-correlation  function  (ACF).  The  process  may  be  deterministic,  stochastic, 
finite-energy,  infinite-energy,  non-time-varying  (stationary)  or  time-varying 
(non-stationary). 

The  ACF  of  a  stochastic  process  is  the  correlation  between  two  samples  of  the 
process  taken  at  ti  and  t2,  and  is  defined  as: 

Rxx(ti,t2)  =  E{x  (t,)  x*(t2)j,  (mi) 

where  E(.)  is  the  expectation  operator  and  *  stands  for  the  complex  conjugation.  For  a 
stationary  (or  wide-sense  stationary)  process,  R(tht2)  depends  only  on  the  time  lag 
r  =  r,  -r2,  resulting  in  a  stationary  spectrum.  The  Wiener-Khinchin  theorem  defines  the 
relationship  between  the  correlation  function  and  spectral  density  as 

5«(®)  =  (T)e~J0*dr.  (m.2) 

B.  INSTANTANEOUS  CORRELATION  FUNCTION  DEFINITION 

The  ACF  of  deterministic  and  stochastic  processes  are  computed  using  time 
domain  averaging  and  the  expectation  operator,  respectively.  This  means  that  a 
smoothing  process  has  to  be  applied  to  compute  the  correlation  functions.  The 
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instantaneous  correlation  function  (ICF)  does  not  use  an  averaging  operation  (integration 
nor  expectation).  The  instantaneous  correlation  function  is  defined  as  the  product  of  two 

samples  of  the  process.  These  two  samples  are  drawn  at  two  time  instants  centered  about 

i 

U  %'rnisimmMiilskMsMdK 

where  i  denotes  the  instantaneous  nature  of  the  correlation  function. 

If  x(t)  is  a  sinusoidal  signal  then  the  multiplication  to  obtain  the  instantaneous  values 
of  R‘(t, T)  generates  cross  terms  in  the  ICF.  For  example,  the  real-valued  sinusoidal 

signal  x(t)=  A  cos  ( cot)  has  an  ACF  given  by: 

A2 

R(  T)  = COS(CUT), 

2 

while  the  ICF  is  given  by: 

R‘  ( t,T )  =  —  [cOS(2t»0  +COS(£0T)]. 

The  ACF  of  a  single  sinusoidal  signal  has  only  one  component  and  no  cross  term, 
while  the  ICF  has  cross  terms.  If  the  signal  x(t)  is  represented  by  its  analytic 

form  x(t)  =  Aejart ,  then  its  ICF  is  given  by. 

R‘(t,T)  =  A2eim. 

That  is,  the  ICF  of  a  single  complex  exponential  signal  has  no  cross-terms. 
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C.  HILBERT  TRANSFORM 

The  Hilbert  transform  allows  the  use  of  the  analytic  signal.  In  general  taking  the 
real  valued  signal  and  adding  the  Hilbert  transform  as  the  imaginary  part  accomplishes 
the  analytic  signal  generation,  i.e.,  one-sided  spectral  density.  We  note  that  in  MATLAB 
taking  the  Hilbert  transform  amounts  to  generating  the  analytic  signal. 

D.  MEDIAN  FILTERING 

Median  filtering  is  a  non-linear  filtering  technique  that  ranks  the  data  in  amplitude 
over  a  window  of  consideration  and  replaces  the  center  of  the  window  with  the  center 
(median)  of  the  ranked  data  vector. 
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IV.  SIGNAL  ANALYSIS  TECHNIQUES  AND  SIMULATION 

A.  PULSE  GENERATION 

For  simulation  purposes  a  frequency  hopped  radar  pulse  sequence  (i.e.,  interpulse 
modulation)  is  generated,  consisting  of  ten  pulses.  At  each  signal  to  noise  (SNR)  level 
100  realizations  are  obtained,  making  a  set  of  1000  pulses  available  at  each  SNR  level. 
The  spacing  between  the  ten  pulses  is  randomized.  The  pulses  can  take  on,  in  a  random 
fashion,  one  of  six  possible  carrier  frequencies.  The  signals  are  used  to  obtain  information 
about  the  false  alarm  rate  (Pfa  ),  the  probability  of  detection  (PD),  and  what  we  define  as 
the  probability  of  multiple  detection  (Pmul)-  The  definition  of  Pfa  and  PD  are  the  classical 
ones,  relating  to  the  probability  of  detecting  a  pulse  when  there  is  noise  only  and  the 
probability  of  detecting  the  pulse  in  the  presence  of  noise,  respectively.  Pmul  denotes  the 
improper  detection  of  more  than  one  pulse  when  only  one  pulse  is  present.  This  can  be  a 
problem  in  an  automated  detection  system  operating  at  low  SNR  levels  or  at  high  SNR 
levels  if  the  threshold  is  set  high.  The  experimental  SNR  ranges  from  3  to  20  dB.  For 
perfect  detection  and  for  perfect  false  alarm  rate  experiments  SNR  values  from  6  to  20  dB 
are  used. 

B.  DETECTION  TECHNIQUES 

Figure  IV.  1  shows  the  four  processing  techniques  described  in  this  chapter.  The 
test  signals  are  real  valued,  but  when  appropriate  their  analytic  signal  versions  are 
processed  [48].  Time  domain  data  in  real  or  analytic  form  is  used  to  obtain  the  time 
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correlation  function  (TCF)  or  it  is  used  directly  in  an  envelope  detector.  The  triangular 
sections  of  the  TCF,  indicating  the  time  ranges  when  pulses  are  present,  have  noise  and 
modulated  signal  components,  while  the  remainder  of  the  surface  has  noise  only.  In  the 
description  of  the  operation,  i.e.,  Figure  IV.  1,  the  block  labeled  “Hilbert”  denotes  the 
operation  to  obtain  the  analytic  signal.  A  detailed  discussion  of  the  TCF  can  be  found  in 
[46,47], 


1.  Real  Processing: 

Real  valued  data  is  used  in  the  computation  of  the  TCF.  To  extract  the 
envelope  of  the  TCF  the  absolute  value  of  its  analytic  signal  representation 
(denoted  by  Hilbert  and  the  summation  symbol)  is  taken.  Next,  the  values  are 
summed  over  both  the  +45°  and  -45°  degree  directions.  The  next  processing  step 
is  a  median  filtering  operation  that  smoothes  the  data  but  retains  the  sharp  edges 
indicating  the  onset  and  the  end  of  a  pulse.  The  final  procedure,  used  as  the  final 
step  in  all  four  processing  techniques,  is  a  threshold  operation.  The  onset  of  a 
pulse  and  end  of  a  pulse  are  declared  when  the  threshold  is  crossed  in  an  upward 
or  downward  fashion,  respectively. 

2.  Triangle  Processing: 

This  operation  is  identical  to  that  described  for  the  Real  Processing  technique 
with  one  exception.  The  exception  is  that  the  Hilbert  transform,  necessary  for 
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obtaining  an  analytic  signal  is  taken  along  the  45°  angle,  rather  than  the  pre¬ 
programmed  90°  directions. 

3.  Complex  Processing : 

For  this  processing  technique  the  analytic  signal  representation  of  the  time 
domain  data  is  used.  From  the  analytic  (complex  valued)  data  the  TCF  is 
computed.  This  operation  is  followed  by  sequentially  taking  the  magnitude 
(absolute  value)  and  the  sum  as  described  in  the  Real  Processing  technique,  i.e., 
summing  up  the  magnitudes  in  both  the  plus  45°  and  minus  45°  degree  directions. 
The  next  step  is  the  median  filtering  operation  to  smooth  the  data.  The  final 
operation  is  the  threshold  operation.  When  the  threshold  is  exceeded  the  onset  of  a 
pulse  is  declared.  Conversely,  the  end  of  a  given  pulse  is  declared  when  falling 
below  the  threshold  level. 

4.  Envelope  Processing : 

As  the  name  implies  the  envelope  is  extracted  by  taking  the  absolute  value  of 
the  analytic  data  (absolute  value  of  the  analytic  signal  values).  The  envelope  of  the 
data  is  filtered  by  a  median  filter,  which  is  followed  by  the  threshold  detection 
mentioned  before. 
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Real  Signal  Reai 


Figure  IV.l:  Envelope-based  and  TCF-based  detectors. 

Figure  IV.  1  shows  the  four  processing  techniques,  outlined  before,  in  a  block  diagram 
form.  The  top  channel  shows  the  “Real  Processing”  scheme,  while  the  second  channel 
shows  the  “Triangle  Processing”  technique.  The  third  and  fourth  channel  show  the 
Complex  and  Envelope  Processing  techniques,  respectively  used. 

C.  SIMULATION  RESULTS 

Threshold  sensitivity  can  be  established  by  examining  the  number  of  pulses 
detected  and  the  threshold.  The  threshold  is  expressed  in  integer  multiples  of  the  noise 
only  standard  deviation.  Figures  IV.2  to  IV.4  show  results  for  SNR  values  of  20dB,  lOdB, 
and  6dB.  In  each  graph  the  minimum  threshold,  in  integer  multiples  of  the  noise  only 
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standard  deviation,  which  allows  a  false  alarm  rate  of  zero  is  provided.  Each  figure 
consists  of  four  graphs,  where  each  graph  shows  simulation  results  for  one  detector. 
There  are  four  detection  schemes,  hence  there  are  four  graphs  per  figure.  The  SNR  is 
fixed  at  a  given  value  for  the  four  graphs  of  a  given  figure. 

Each  graph  can  be  broken  into  four  sections.  In  the  first  section  the  threshold  is 
too  low,  hence  the  noise  causes  false  alarms.  If  the  threshold  is  increased  above  the  noise 
peaks  only  the  signal  is  detected.  If  the  threshold  is  further  increased  a  region  where  the 
processed  pulses  begin  to  break  apart,  is  reached.  This  causes  the  number  of  detected 
pulses  to  increase  (a  given  pulse  becomes  multiple  events,  or  multiple  pulses).  In  the 
fourth  region  the  number  of  detection  begins  to  drop,  eventually  going  to  zero  as  the 
threshold  is  increased. 

Each  graph  displays  two  curves.  The  lower  one  shows  the  probability  of  detection, 
while  the  upper  one  shows  the  total  number  of  detection,  expressed  in  a  normalized 
percentage.  For  example,  a  value  of  2.2  for  the  top  trace,  at  a  value  of  0.9  for  the  lower 
trace,  corresponds  to  a  total  number  of  detection  of  220  percent  with  a  probability  of 
detection  of  90%  and  a  number  of  (improper)  multiple  detection  of  130%.  So  for  1000 
realizations,  one  has  900  proper  detections,  1300  improper  multiple  detections  and  hence 
a  total  of  2,200  detections.  This  information  is  useful  in  assessing  an  appropriate 
threshold  setting  and  to  determine  over  what  range  of  thresholds  one  can  accurately  detect 
the  number  of  radar  pulses. 
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Figure  2a.  Performance  curve  normalized  percentage  versus  noise  sigma,  SNR=20dB. 

Figure  IV.2:  Performance  curve;  normalized  percentage  versus  noise  sigma;  SNR=20dB. 
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Figure  2b.  Performance  curve  normalized  percentage  versus  noise  sigma,  SNR=10dB. 

Figure  IV.3:  Performance  curve;  normalized  percentage  versus  noise  sigma;  SNR=10dB 


Figure  2c.  Performance  curve  normalized  percentage  versus  noise  sigma.  SNR=6dB. 

Figure  IV.4:  Performance  curve;  normalized  percentage  versus  noise  sigma;  SNR=6dB. 
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Table  IV.  1  uses  a  threshold  of  three  noise  standard  deviations  for  each  of  the  four 
detectors.  It  shows  the  performance  in  terms  of  the  Pfa  as  a  function  of  SNR  for  a  Pd  of 
1.0,  except  when  specially  denoted.  It  is  apparent  that  the  Real  (followed  by  the 
Complex)  Processing  scheme  is  the  most  promising  one  in  the  SNR  range  between  6  and 
20  dB. 


Pd  =  1.0 

SNR 

20  dB 

16  dB 

13  dB 

10  dB 

6  dB 

3  dB 

TRIANGLE 

0.046 

0.056 

0.060 

0.095 

0.165 

0.236 

COMPLEX 

o" 

0 

0.001 

0.012 

0.047 

0.120 

REAL 

0 

0 

0 

0.004 

0.039 

0.091  @ 
Pfa=0.99 

5 

ENVELOPE 

0 

0.006 

0.022 

0.045 

0.094 

0.133  @ 
Pfa=0.99 

8 

Table  IV.l:  Performance  for  fixed  threshold. 

Table  IV.2  shows  the  performance  in  terms  of  the  Pd  as  a  function  of  SNR  for  a 
PFA  equal  to  0.  It  is  apparent  that  the  Complex  (followed  by  the  Envelope)  Processing 
scheme  is  the  most  promising  one  in  the  SNR  range  between  3  to  20  dB. 
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Pfa  = 

0.0 

SNR 

20  dB 

16  dB 

13  dB 

10  dB 

6  dB 

3  dB 

0  dB 

TRIANGLE 

1.0 

1.0 

1.0 

1.0 

0.973 

0.438 

0.036 

COMPLEX 

1.0 

1.0 

1.0 

1.0 

1.0 

0.653 

0.029 

REAL 

1.0 

1.0 

1.6 

i.o 

0.9677 

0.479 

0.071 

ENVELOPE 

1.0 

1.0 

1.0 

1.0 

1.0 

0.312 

0.017 

Table  IV.2:  Performance  for  fixed  PFa  =  0. 

Table  IY.3  allows  the  comparison  of  the  detectors  in  terms  of  the  sensitivity  to  the 
multiple  detection  of  a  given  pulse.  For  the  tests  executed  for  this  section,  100 
realizations  of  a  pulse  train  with  10  members  is  undertaken.  This  fixes  the  total  number  of 
possible  pulses  to  1000  for  each  SNR.  The  top  line  shows  the  number  of  properly 
detected  pulses,  while  the  lower  line  shows  the  number  of  multiple  detection  of  the  pulse 
trains.  It  is  apparent  that  the  Real,  followed  by  the  Envelope,  detectors  provide  the  most 
robust  schemes. 
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Pd=  10 


20  dB  16  dB  13  dB 

1000  1000  1000 

0 _ 0 _ 0_ 

1000  1000  1000 

0 _ 0 _ 0_ 

1000  1000  1000 

0 _ 0 _ 0_ 

1000  1000  1000 

o  1  0  I  0 

Table  IV.3:  Detection  versus  SNR  level. 

It  becomes  apparent  that  with  pulsed  signals,  depending  on  the  importance  of  the 
undesirable  characteristics,  any  of  the  three  detectors  (Complex,  Real  or  Envelope)  can  be 
used.  From  the  point  of  view  of  processing  cost,  the  Envelope  detector  is  preferable  over 
the  other  two.  From  the  point  of  false  alarm  rate  or  probability  of  detection  the  Complex 
or  Real  detector,  respectively  is  more  desirable. 


SNR 

TRIANGLE 

COMPLEX 

REAL 

ENVELOPE 


10  dB 
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V.  TIME-FREQUENCY  REPRESENTATIONS 


A.  BACKGROUND 

Numerous  schemes  have  been  proposed  over  the  years  to  retrieve  modulation 
parameters  of  chirp  signals  for  various  types  of  applications  ranging  from 
communications  to  radar  and  sonar.  This  chapter  focuses  on  signals  with  constant 
amplitudes  and  linear  and  hyperbolic  intra-pulse  modulation.  Such  types  can  be  found  in 
pulse-compression  radar  signals,  pulse  or  CW  radar  signals  received  from  targets,  or 
modulation  signals  found  in  communication  applications.  Radar  intra-pulse  modulation 
parameters  are  usually  extracted  using  hardware  schemes,  as  they  are  well  suited  to 
extract  the  instantaneous  frequency  at  high  SNR  levels.  In  this  work  we  consider  the 
problem  from  a  different  angle  and  investigate  the  application  of  time-frequency  (TF)  and 
basic  image  processing  techniques  to  extract  the  modulation  information. 

This  work  assumes  that  1)  individual  pulses  have  been  detected  and  isolated  prior 
to  the  processing  discussed  next,  and  2)  the  modulation  type  is  known,  as  we  do  not 
discuss  issues  related  to  the  specific  identification  of  the  modulation  scheme  which  is  left 
for  further  research. 

B.  INTRODUCTION 

Several  different  types  of  approaches  to  estimate  the  modulation  parameters  from 
chirp  signals  have  been  proposed  in  the  literature  over  the  years  [25-31,33,34,36-40]. 
Maximum  likelihood  estimates  of  the  chirp  parameters  can  be  derived  in  theory. 


23 


However  they  often  lead  to  complicated  algorithms,  which  may  require  n-dimensional 
optimization  steps  depending  on  the  specific  modulation  type.  Phase  unwrapping 
techniques,  application  of  the  high-order  ambiguity  function  [27],  sub-optimal  schemes  to 
the  maximum  likelihood  estimates  [25],  and  others  have  also  been  proposed,  each  with  its 
own  advantages  and  drawbacks.  In  addition,  combinations  of  TF  representations  and 
pattern  recognition  schemes  have  been  considered  to  estimate  chirp  parameters  [28- 
30,33,36]. 

Our  study  considers  the  combination  of  TF  representations  and  pattern  recognition 
schemes  to  estimate  the  chirp  modulation  parameters,  and  compare  the  resulting 
performances  in  various  white  Gaussian  noise  levels.  The  schemes  proposed  by 
Barbarossa  et  al.  [28,29]  and  Wood  and  Barry  [36]  both  combine  the  TFR-based  Wigner- 
Ville  (WV)  transformation  and  the  Radon  (Hough)  transform  to  estimate  the  chirp 
parameters,  as  the  WV  transform  is  well  localized  for  linear  chirps.  Barbarossa  also 
indicated  that  any  TFR  could  potentially  be  considered,  provided  that  it  has  good 
localization  and  concentration  properties,  and  is  robust  to  noise  degradations.  However, 
no  performance  comparisons  were  made  available. 

Our  study  considers  a  wide  range  of  TFRs  and  identifies  the  top  three 
representations  best  suited  to  the  follow-on  extraction  scheme.  We  selected  eleven 
different  time-frequency  and  wavelet  transformations  (TFR)  and  applied  a  Radon-based 
transform  step  to  extract  the  parameters  from  linear  chirps,  as  discussed  in  Sections  VI.B 
and  VI.C.  The  selection  of  the  top  transformations  was  made  by  comparing  the  resulting 
errors  in  the  estimation  of  the  modulation  parameters,  as  discussed  in  Section  VI.C.3. 
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The  second  phase  of  the  study  considered  the  extraction  of  modulation  parameters 
for  hyperbolic  chirps.  We  restricted  our  study  to  the  top  three  energy  distributions  leading 
to  the  best  image  quality  for  the  linear  chirp,  and  derived  the  estimation  procedure 
presented  in  Section  VII. 

C.  TIME-FREQUENCY  REPRESENTATIONS  DESCRIPTION 
1.  Introduction 

Generally,  a  signal  can  be  represented  using  different  types  of  decompositions. 

The  two  most  widely  used  representations  are  the  time  domain  and  the  frequency  domain 

representation.  The  first  shows  how  the  amplitude  of  the  signal  changes  with  respect  to 

time  while  the  second  shows  how  often  these  changes  occur.  However,  the  Fourier 

transform  is  not  well  suited  to  represent  non-stationary  signals,  as  it  provides  no  exact 

information  regarding  the  variations  of  the  signal  characteristics  as  a  function  of  time. 

Analyzing  such  signals  requires  a  joint  time-frequency  representation.  Several  techniques 

have  been  developed  for  this  type  of  representation.  This  study  considered  the  following 

eleven  transformations  from  the  two  major  classes  of  atomic  decompositions  and  energy 

distributions.  The  goal  was  to  select  a  small  number  of  transformations  leading  to  the  best 

“image  quality”  from  that  set.  We  briefly  introduce  each  of  them.  Further  details 

regarding  each  of  them  may  be  found  in  Moraitakis  [33,  Sect.  HI].  The  transformations 

considered  are: 

-Wavelet  packet  best  basis, 

-Cosine  packet  best  basis, 

-Wavelet  pursuit, 

25 


-Cosine  pursuit, 

-Wigner-Ville  distribution, 

-Pseudo  Wigner-Ville  distribution, 

-Reassigned  Pseudo  Wigner-Ville  distribution, 

-Smoothed  Pseudo  Wigner-Ville  distribution, 

-Reassigned  Smoothed  Pseudo  Wigner-Ville  distribution, 

-Spectrogram  distribution, 

-Reassigned  spectrogram  distribution. 

For  comparison  purposes,  Figure  V.l  shows  the  resulting  images  obtained  for  a  linear 
chirp  with  SNR=10dB  when  no  energy  thresholding  is  applied  to  the  images. 
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Figure  V.l:  Various  time-frequency  representations,  linear  chirp,  SNR=10dB,  Frequency 
normalized  to  sampling  frequency. 
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2.  Description 

Wavelet  and  cosine  decompositions  were  implemented  with  the  software 
“Wavelab  7.01”  software  [12],  while  all  others  used  the  “Time-Frequency  Toolbox” 
[9,11]. 

1 .  Wavelet  packet  best-basis  decomposition  (WPD ): 

This  decomposition  performs  a  multi-resolution  decomposition  by  partitioning  the 
frequency  axis.  It  can  be  viewed  as  an  extension  of  the  wavelet  transform  where  both 
high-pass  and  low-pass  components  are  decomposed,  thereby  allowing  for  more 
flexibility  in  the  decomposition  than  the  wavelet  transform  does.  Selection  of  the  specific 
basis  (called  the  “Best  Basis”)  is  obtained  by  using  a  user-specified  information  cost 
criterion  [3,5-7]. 

2.  Cosine  packet  best-basis  decomposition  (CPD): 

This  decomposition  is  similar  in  concept  to  the  WP  decomposition,  but  performs  a 
multi-resolution  decomposition  by  partitioning  the  time  axis  instead.  It  is  well  matched 
to  narrowband  signals  [5]. 

3.  Wavelet  pursuit  (WP): 

The  wavelet  pursuit  is  an  iterative  procedure  where  the  signal  is  decomposed 
adaptively  into  a  set  of  functions,  not  necessarily  orthonormal.  This  procedure  has  been 
shown  to  be  well  suited  to  represent  non-stationary  signals  with  varying  time-frequency 
characteristics  [4,8], 
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4.  Cosine  pursuit  ( CP): 

The  cosine  pursuit  is  similar  in  concept  to  the  wavelet  pursuit,  where  the  signal  is 
decomposed  into  a  set  of  locally  defined  cosine  functions  [4]. 


5.  Wigner-Ville  distribution  (WVD): 

The  WVD  belongs  to  the  Cohen’s  class  of  energy  distributions  [2,10].  It  allows 
for  perfect  localization  of  linear  chirps.  However,  its  quadratic  definition  is  responsible 
for  interference  terms  that  appear  when  the  signal  is  not  linearly  modulated  or  more  than 
one  signal  is  present.  Figure  V.2  shows  the  WVD  and  the  effect  of  the  cross  terms  on  four 
different  signals.  The  effect  is  negligible  in  the  noise-free  linear  chirp  scenario  only. 
Cross-terms  show  as  a  line  between  the  two  true  components  when  the  signal  is 
composed  of  two  parallel  chirps. 


Figure  V.2:  Wigner-Ville  distribution;  time  normalized  to  the  pulse  length. 
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6.  Pseudo  Wigner-Ville  distribution  (PWVD): 

The  PWVD  introduces  a  time-window  to  the  WVD  definition,  thereby  reducing 
the  interferences  present  in  the  transform,  as  illustrated  in  Figure  V.3.  However,  the  time 
windowing  worsens  the  frequency  resolution  and  the  PWVD  looses  many  of  the  valuable 
properties  of  the  WVD  [2]. 


Two  linear  Chirps 


Linear  Chirp  in  Noise  (SNR=2db) 


Figure  V.3:  Pseudo  Wigner-Ville  distribution;  time  is  normalized  over  the  pulse  duration. 


7.  Smoothed  Pseudo  Wigner-Ville  distribution  (SPWVD): 

The  SPWD  provides  an  attenuation  of  the  interference  terms  by  smoothing  both  in 
in  frequency  and  time  domains.  The  drawback  of  this  method  is  that  the  frequency 
resolution  further  degrades.  Figure  V.4  shows  that  the  cross-terms  have  almost  been 
eliminated  while  the  frequency  resolution  has  degraded. 
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Two  linear  Chirps 


Hyperbolic  Chirp 


Figure  V.4:  Smoothed  Pseudo  Wigner-Ville  distribution;  time  is  normalized  over  the  pulse 
duration. 

8.  Spectogram: 

The  spectrogram  is  defined  as  the  magnitude  of  the  Short-Time  Fourier  transform 

[2]. 


9.  Reassigned  spectogram: 

We  see  that  the  WVD  resolution  worsens  when  we  try  to  suppress  the  interference 
terms.  The  reassignment  method  was  introduced  in  an  attempt  to  improve  the  resolution 
of  the  transformation  [9].  This  scheme  assumes  that  the  energy  distribution  in  the  time- 
frequency  plane  resembles  a  mass  distribution  and  moves  each  value  of  the  time- 
frequency  plane  located  at  a  point  (t,f)  to  another  point  (t’,f ),  which  is  the  center  of 
gravity  of  the  energy  distribution  in  the  area  of  (t,f).  The  result  is  a  very  focused 
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representation  with  high  intensity  since  the  value  at  the  point  (t’,f )  is  the  sum  of  all  the 
neighboring  values.  The  reassignment  method  can  be  applied  to  most  energy 
distributions  as  well  as  to  the  spectrogram.  However,  the  computational  load  is  quite 
high.  The  top  plot  in  Figure  V.5  shows  the  hyperbolic  chirp  used  in  Figures  V.2  to  V.4,  as 
represented  for  the  reassigned  spectrogram  distribution. 


Reassigned  Spectogram 


0  .  1 


Figure  V.5:  Reassignment  method  applied  to  a  hyperbolic  chirp  test  signal;  time  is 
normalized  over  the  pulse  duration. 

10.  Reassigned  Pseudo  Wigner-Ville  distribution  (RPWVD): 

The  reassignment  method,  as  discussed  above,  is  applied  to  the  PWVD.  The 

middle  plot  in  Figure  V.2  shows  the  hyperbolic  chirp  used  in  Figures  V.2  to  V.4,  as 

represented  for  the  reassigned  RPWVD. 
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11.  Reassigned  Smoothed  Pseudo  Wigner-Ville  distribution  (RSPWV): 

The  reassignment  method,  as  discussed  above,  is  applied  to  the  SPWVD.  The 

bottom  plot  in  Figure  V.5  shows  the  hyperbolic  chirp  used  in  Figures  V.2  to  V.4,  as 
represented  for  the  reassigned  RSPWVD.  Note  the  quality  of  the  focusing  in  all  three 
reassigned  methods  and  the  resulting  lack  of  cross  terms. 

Simulations  showed  that  the  best  image  quality  for  the  various  time-frequency  and 
time-scale  representations  of  the  noisy  chirps  is  obtained  with  the  time-frequency 
transformations,  as  illustrated  in  Figure  V.l  above  which  shows  a  more  focused 
representation  of  the  information.  Therefore,  one  would  expect  better  estimation  of  the 
chirp  parameters  when  using  time-frequency  representations  than  with  time-scale 
represenations  as  the  SNR  decreases.  Section  VI  confirms  these  findings  by  investigating 
the  estimation  performance  obtained  from  each  of  the  time-frequency/scale 
transformations. 
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VI.  APPLICATION  OF  THE  RADON  TRANSFORM  TO  LINEAR 
CHIRP  PARAMETERS  ESTIMATION 

Once  the  time-frequency  representation  of  the  signal  information  is  obtained,  the 
Radon  transform  is  applied  to  extract  the  modulation  parameters  for  linear  chirp  signals. 
The  Radon  transform  has  been  used  extensively  in  image  processing  for  edge  detection 
[13,14]  and  has  been  generalized  to  detect  curves  of  arbitrary  shapes  [17].  We  briefly 
review  its  concept,  before  presenting  its  application  to  our  problem. 

A.  INTRODUCTION 

Let’s  assume  we  have  an  arbitrary  function  f(x,y)  in  a  subspace  of  R2.  The  two- 
dimensional  Radon  transform  is  defined  as  the  projection  or  line  integral  of  the  function 
f(x,y)  along  all  possible  lines  L  [13].  Mathematically,  the  transformation  is  described  as: 

R  =  jf(x(s,L),y(s,L))ds.  (VI-1) 

L 

Recall  that  the  equation  of  a  line  in  polar  coordinates  is  given  by: 

p  =  x  •  cos($)  +  y  •  sin($) ,  (VI-2) 

where  p  and  $  represent  the  distance  from  the  origin  and  the  angle  measured 
counterclockwise  from  the  x  axis  respectively,  as  shown  in  Figure  VI.  1.  Now  suppose  we 
use  another  coordinate  system  with  axes  rotated  by  the  angle  #.  The  new  x-axis  lies  on 
the  line  with  associated  orthogonal  direction  “s”.  The  two  cartesian  systems  xoy  and  pos 
are  related  to  each  other  via  the  following  relation: 
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X 


(VI-3) 


Equation  (IV- 1)  can  be  rewritten  as: 


i?(p,tf)=  J  f  (p  cost?  -  s  sin  t),  p  sin  t)  +  s  cos  t))  •  ds 


(VI-4) 


The  above  equation  shows  that  the  Radon  transform  translates  a  two-dimensional 
function  of  the  variables  (x,y)  to  one  with  variables  (p,t)) .  Thus,  the  Radon  transform  of 
an  image  taken  at  a  specific  angle  $  is  the  projection  of  the  image  onto  the  line  which 
forms  an  angle  t)  with  the  x-axis. 


Figure  VI.l:  Two-dimensional  Radon  transform. 

The  Radon  transform  of  a  single  line  with  a  slope  angle  (p  for  the  specific  angle 
7T 

$  =  — t-  cp  is  a  single  point  with  intensity  equal  to  the  sum  of  the  intensities  at  each  point 
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on  the  line.  This  property  allows  detection  of  lines  in  an  image.  In  addition,  the  transform 
is  also  robust  to  noise  degradations. 


B.  LINE  PARAMETER  IDENTIFICATION 

Line  parameters  can  be  obtained  using  the  Radon  transform  as  follows:  Assume 

the  image  under  investigation  contains  a  line  which  has  an  angle  q>  with  the  x-axis,  as 
shown  in  Figure  VI.2.  The  equation  of  the  line  can  be  defined  in  terms  of  its  slope  a  and 
initial  offset  value  b. 

y  =  a-x  +  b,  a  =  tan(<p) .  (VI-5) 

The  Radon  transform  R(p,i3)  of  the  image  for  angles  0°  to  180°  is  maximum  when 
the  projection  of  the  line  has  a  minimum  area.  Thus,  it  is  maximum  when  the  line  is 
perpendicular  to  the  projection  line,  i.e.,  when  $=(p+90",  which  leads  to 


a  =  tan(tf-90°). 


(VI-6) 


Now,  if  we  take  the  Radon  transform  of  the  image  for  the  specific  angle 
t9=<p+90°,  we  can  estimate  the  offset  parameter  b  from  the  position  of  its  maximum  value 


along  the  axis  p.  This  position  is  indicated  in  Figure  VI.2  as  “C”.  Next,  the  offset 
parameter  b  can  be  computed  as: 


OC-^s  in(tf-90“)  „ 

fc  =  — + - £ -  =  — -  +  - 

2  cos(#-90°)  2 


OC  +  — -^-cos(?7) 
2 

sin($) 


(VI-7) 
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where  Ny  is  the  length  of  the  vertical  side  of  the  image,  while  N  is  the  length  of 
the  horizontal  side  (in  this  case  Nx=Ny=256).  Note  that  the  distance  OC  can  be  positive  or 
negative  and  is  actually  negative  in  Figure  VI.2. 

Thus,  the  equation  of  the  line  can  be  determined  by  first  applying  the  Radon 
transform  of  the  image  for  angles  between  0°  to  180°,  then  finding  the  coordinates  of  the 
maximum  point  (tSU,,  Pmax),  as  shown  in  Figure  VI.3,  and  finally  using  equations  (VI-6) 
and  (VI-7)  using  OC=pmax  and  t?=tW  Note  that  the  accuracy  of  the  estimates  depends 
on  the  resolution  of  the  image  and  the  size  of  the  angle  step  selected  for  the  Radon 
transform.  Unfortunately,  any  attempt  to  increase  the  resolution  (size)  or  decrease  the 
angle  step  results  in  increasing  the  computation  load. 


Figure  VI.2:  Radon  transform  geometry;  image  containing  a  single  line. 
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Figure  VI.3:  Radon  transform  for  line  shown  in  Figure  VI.2;  for  $=0°,...,180°,  1°  inc. 

To  minimize  the  computational  load  we  apply  the  Radon  transform  in  two  stages. 
First,  we  scan  the  image  from  0°  to  179°  in  steps  of  2°.  We  determine  the  angle  -dim  that 
corresponds  to  the  column  that  includes  the  point  with  the  highest  intensity.  Next,  we 
scan  the  image  from  the  angle  fhm-l°  to  fhm+l°  in  steps  of  0.1°.  Using  this 
implementation  allows  the  maximum  line  slope  error  to  decrease  to  ±0.05°  without 
having  to  cover  the  whole  range  of  angles  at  that  rate. 

C.  LINEARLY  MODULATED  CHIRPS 

This  section  presents  the  scheme  used  to  estimate  linear  chirp  parameters.  First, 
we  present  the  method  used  to  generate  the  signals.  Next,  we  describe  the  image 
processing  technique  applied  to  extract  the  features  from  the  time-frequency  image. 
Finally,  we  compare  the  results  obtained  for  each  time-frequency  representation 
considered  and  their  robustness  to  noise  degradations. 
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1.  Signal  Generation 

The  radar  signals  considered  in  this  study  are  synthetic  and  consist  of  a  train  of 
several  linearly  modulated  pulses,  as  shown  in  Figure  VI.4. 


Figure  VI.4:  Synthetic  radar  signal;  linear  chirp  modulation. 


First,  we  assume  that  we  can  isolate  one  single  pulse  of  duration  x.  This  extraction 
can  easily  be  done  with  an  energy  detector  in  medium  to  high  SNR  levels.  Thus,  the 
received  signal  has  an  instantaneous  frequency  f(t)  defined  as: 


f(t)  =  fo+k-t  (VI-8) 

Some  of  the  time-frequency  transformations  considered  in  the  study  use  the 
analytical  version  derived  from  the  real  received  signal.  In  such  cases,  the  analytical 
signal  is  computed  with  the  Hilbert  transform.  The  analytical  noise-free  linear  chirp 
signal  s(t)  is  given  as: 


s(t)  =  e 


(VI-9) 
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The  signal  is  assumed  to  be  sampled  a  rate  of  512  [samples/pulse  length].  Recall  that  the 
sampling  frequency  for  a  real  signal  must  be  equal  to  at  least  twice  its  highest  frequency 
component.  Thus,  the  signal  may  be  heterodyned  down  to  a  lower  frequency  if  needed. 
The  discrete  real  part  of  the  signal  is  given  as: 

J2x(^.n+~n2)  f  Tr 

5[n]  =  Re(e  u  2f‘  )  =  cos(2  n(^-n  +  — T-n2)),  n=l,..,512  (VI-10) 

/,  ~  ft 

We  can  rewrite  the  above  equation  as: 

s[n]  =  cos(2tt(/m0 -n  +  ^-n2))  where  fnQ  =  ^-,  a  =  ,  n=l,..,512 

2  fs  fs 

(VI-11) 

The  terms  a  and  fno  represent  the  normalized  slope  and  the  normalized  starting 
frequency  respectively,  with  respect  to  the  sampling  frequency  fs.  Thus,  the  normalized 
frequency  equation  for  the  linear  chirp  discrete  signal  is  given  by: 

fn  =  /«o  +a-n  ,  n=l,..,512.  (VI-12) 

Normalizing  the  time  index  over  the  pulse  length  leads  to  the  following 
normalized  frequency  equation  for  the  linear  chirp: 

fnn=fn0/N  +  (a/N)-n,  n=l,..,512. 

Multiple  linear  chirp  signals  trials  were  generated  by  randomly  selecting  both  the 
initial  frequency  and  the  slope.  Note  that  sampling  constraints  need  to  be  satisfied  to 
avoid  aliasing  in  the  resulting  discrete  linear  chirp.  For  example,  the  parameter  a  needs  to 
be  selected  so  that  the  final  instantaneous  frequency  doesn’t  exceed  0.5  for  a  given  initial 
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frequency /„o.  As  a  result,  aliasing  may  be  avoided  by  selecting  both  the  initial  and  final 
normalized  instantaneous  frequencies  within  the  range  [0,  0.5],  and  computing  the 
corresponding  slope  parameter  a.  Initial  and  final  frequencies  are  selected  from  a 
uniform  distribution  defined  in  the  range  [0,0.5].  Next,  the  chirp  is  corrupted  by 
additive  white  Gaussian  noise.  Analytical  expressions  needed  to  compute  some  of  the 
time-frequency  transformations  are  obtained  with  the  Hilbert  transform  of  the  noisy 
signal.  Real  signal  expressions  were  used  to  compute  the  wavelet-based  decompositions. 


2.  Simulation  Set-up  and  Extraction 

The  eleven  different  time-frequency  and  wavelet-based  representations  described 

in  Section  V  were  considered.  The  goal  was  to  select  a  small  number  of  transformations 

leading  to  the  best  “image  quality”  from  that  set.  Recall  that  the  representations 

considered  were: 

-Wavelet  packet  best  basis, 

-Cosine  packet  best  basis, 

-Wavelet  pursuit, 

-Cosine  pursuit, 

-Wigner-Ville  distribution, 

-Pseudo  Wigner-Ville  distribution  (PWVD), 

-Reassigned  Pseudo  Wigner-Ville  distribution  (Reassigned  PWVD), 

-Smoothed  Pseudo  Wigner-Ville  distribution  (SPWVD), 

-Reassigned  Smoothed  Pseudo  Wigner-Ville  distribution  (Reassigned  SPWVD), 
-Spectrogram, 

-Reassigned  spectrogram. 
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Figure  VI.5  presents  the  image  obtained  for  the  various  time-frequency 
representations  considered  for  a  linear  chirp  with  SNR=10dB.  A  few  comments  are  in 
order  regarding  the  selection  of  the  transform  parameters. 

1.  Number  of  Atoms 

Recall  that  defining  a  line  in  a  plane  requires  two  points  only.  However,  these  two 
points  must  be  perfectly  localized  in  both  time  and  frequency,  and  must  be  immune  from 
any  noise  interference.  Theoretically,  we  could  use  two  atoms  only  from  an  atomic 
decomposition  to  define  the  linear  frequency  equation.  However,  this  doesn’t  hold  in 
practice,  as  the  atoms  are  not  perfectly  localized.  Note  that  a  larger  number  of  atoms  may 
better  represent  the  line  trend.  However,  some  of  the  atoms  may  represent  noise 
contributions  for  noisy  signals.  Therefore,  selecting  the  number  of  atoms  to  represent  a 
linear  chirp  in  noise  requires  a  trade-off  between  these  two  issues:  fewer  atoms  to  denoise 
the  signal  and  more  atoms  to  improve  the  resolution.  We  selected  10  atoms  per 
decomposition  for  the  atomic  decompositions  used  in  our  study  after  running  several  trial 
cases. 

2.  Maximum  Decomposition  Level 

Next,  the  maximum  decomposition  level  was  set  to  7,  as  simulations  showed  no 
advantage  in  going  to  higher  levels. 

3.  Wavelet  Type 

The  mother  wavelet  function  was  selected  after  several  trials  among  the  readily 
provided  functions  in  the  Wavelab  software  [12].  Finally,  we  selected  Daubechie-IV 
since  it  gave  the  best  time-frequency  representation  for  the  types  of  signal  considered. 
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4.  Image  Thresholding 

No  image  intensity  thresholding  was  applied  to  the  time-frequency  images,  as 
simulations  showed  that  it  this  step  worsened  the  results  for  low  SNR  levels.  Note  that 
implicit  denoising  is  actually  performed  by  selecting  a  small  number  of  atoms  for  the 
atomic  decompositions. 

5.  Window  Length 

The  PWV,  SPWV  and  Reassigned  SPWVD  transformations  use  a  frequency 
window  which  has  a  Hamming  (N/4)  time  domain  expression,  and  a  Hamming  (N/10)  for 
time  smoothing  window.  The  analysis  window  selected  for  the  spectrogram  is  Hamming 
(N/4),  where  N  is  the  length  of  the  signal  (512  bins). 
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Reassigned  Spectrogram 


Figure  VI.5:  Various  time-frequency  representations  for  a  linear  chirp;  SNR=10dB, 
Normalized  frequency.  Frequency  normalized  to  sampling  frequency. 

6.  Radon  Transformation  Implementation  Issues 

The  Radon  transform  is  selected  to  extract  the  line  equation  from  the  time- 
frequency  image,  as  described  earlier  in  Section  VLB.  We  use  a  two-stage 
implementation  with  final  degree  increment  0.1°.  The  size  of  the  resulting  image  for  each 
time  frequency  representation  is  set  to  256x256  points,  as  the  Radon  transform  for  larger 
images  would  be  too  computationally  expensive  for  a  MATLAB  implementation.  One 
hundred  randomly  generated  linear  chirps  for  a  given  SNR  level  were  generated,  and  the 
SNR  varied  in  the  range  -lOdB  to  lOdB.  Next,  the  chirp  parameters  were  estimated  from 
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the  images  generated  by  each  of  the  eleven  time-frequency  transformations  considered. 
Performance  comparisons  are  presented  next. 


3.  Performance  Results 

3. a)  Evaluation  of  Time-frequency  Representations 

Recall  that  the  maximum  theoretical  slope  error  is  ±0.05°  since  the  final  step 
angle  in  the  radon  transform  is  0.1°,  as  discussed  in  Section  V.  However,  we  also  have  to 
add  noise  effects,  and  quantization  errors  introduced  by  the  finite  resolution  in  the  image. 
Figures  VI.7  to  VI.9  present  the  mean  and  the  standard  deviation  of  the  absolute  slope 
and  offset  errors  as  a  function  of  the  SNR  level  for  all  time-frequency  transformations 
considered  in  the  study.  Figure  VI.6  represents  an  idealized  image  which  could  be 
obtained  with  a  time-frequency  transformation.  Note  that  the  frequency  axis  is 
normalized  by  the  sampling  frequency,  and  the  time  axis  is  normalized  by  the  pulse 
length.  Thus,  the  chirp  slope  parameter  expressed  in  degrees  is  given  by: 

a  = - tan  (/,-/„), 

it 

and  the  offset  value  expressed  in  number  of  (frequency)  bins  is  given  by: 
fnO  =  /oX512, 

where  f0  and  f  respectively  represent  the  normalized  starting  and  ending  chirp 
frequencies.  The  offset  value  is  multiplied  by  the  length  of  the  frequency  transformation 
used. 
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Figure  VI.7:  Slope  and  offset  errors  for  wavelet  and  cosine  packet,  wavelet  and  cosine 
pursuit  decompositions  for  linear  chirp  signals;  SNR  given  in  dB. 
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Figure  VI.9:  Slope  and  offset  errors  for  the  Wigner-Ville  decomposition,  Spectrogram,  and 
reassigned  Spectrogram;  linear  chirp  signals;  SNR  given  in  dB. 

A  few  comments  are  in  order: 

1 .  Results  show  that  all  the  energy  distributions  perform  better  than  the  atomic 
(i.e.,  wavelet-based)  decompositions.  This  is  to  be  expected  as  they  provide  a 
more  accurate  “image”  in  the  time-frequency  plane.  Note  that  atoms  cannot 
be  well  localized  in  both  time  and  frequency,  as  would  be  required  to 
represent  linear  chirps  accurately.  The  best-basis  cosine  packet  decomposition 
gives  the  best  results  followed  by  the  cosine  pursuit  scheme  for  atomic 
decompositions. 
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2.  Most  of  the  energy  distributions  have  slope  errors  very  close  to  the  theoretical 
value  of  0.05°  for  medium  to  high  SNR  levels.  All  transformations  show  a 
SNR  level  breaking  point  at  which  the  errors  suddenly  increase.  This  is  to  be 
expected  as  the  extraction  relies  on  the  quality  of  the  TFR.  Figures  VI.  10  and 
VI.  1 1  illustrate  the  degradation  which  may  occur  in  the  TFR  image  quality  at 
SNR=-10dB  and  the  corresponding  parameter  estimation  results. 

3.  The  Wigner-Ville  distribution  has  a  very  good  and  almost  stable  performance 
for  SNR’s  in  the  range  of  -6dB  to  lOdB.  The  smoothing  time  window  present 
in  the  Pseudo  Wigner-Ville  distribution  improves  the  estimation  in  higher 
SNR  but  shrinks  the  effective  range.  The  Smoothed  Pseudo  Wigner-Ville 
distribution  has  a  slightly  worst  performance  at  high  SNR  but  also  has  the 
widest  effective  range.  The  smoothing  in  time  and  in  frequency  eliminates  the 
interference  terms  almost  completely  so  that  the  representation  is  more 
immune  to  the  noise  than  other  transformations  are.  However,  frequency 
smoothing  results  in  lowered  performance  at  high  SNR  levels. 

4.  Results  show  that  the  reassign  method  usually  improves  the  performance  at 
high  SNR  levels,  as  it  forces  the  representations  to  be  more  “focused”. 
Unfortunately,  applying  the  reassigned  method  in  low  SNR  levels  worsens  the 
performances.  This  is  to  be  expected  as  the  presence  of  noise  close  to  the  line 
moves  the  local  center  of  gravity  of  the  distribution  away  from  its  theoretical 
value. 
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Figure  VI.IO:  Linear  chirp  trial  case  1;  SNR=-10dB;  true  and  estimated  modulation 
parameters. 
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Figure  VI.11:  Linear  chirp  trial  case  2;  SNR=-10dB;  true  and  estimated  modulation 
parameters. 


Table  VI.  1  lists  the  mean  absolute  error  for  the  frequency  slope  and  offset 
parameters  at  SNR  equal  to  lOdB  for  all  transformations  considered  for  a  “rough” 
comparison  of  the  transformation  performances.  In  addition,  Table  VI.  1  lists  the  SNR 
level  at  which  the  errors  suddenly  increase,  indicating  the  SNR  level  at  which  using  a 
specific  distribution  may  become  more  questionable. 


Mean  of  absolute  slope 
error  @  SNR=10dB 
(in  degrees) 

Mean  of  absolute  offset 
error  @  SNR=10dB 
(in  bins) 

Breaking  point 
(in  dB) 

Wavelet  Packet 

0.375 

4 

0 

Cosine  Packet 

0.5149 

3.98 

0 

Wavelet  Pursuit 

0.37 

2.35 

4 

Cosine.  Pursuit 

0.4297 

2.35 

0 

WV 

0.1198 

1.08 

-7 

PWV 

0.0623 

1.02 

-4 

Reassigned  PWV 

0.0729 

1 

-5 

SPWV 

0.0792 

1.02 

-7 

Reassigned 

SPWV 

0.0563 

1.02 

-7 

Spectrogram 

0.098 

1.77 

-7 

Reassigned 

Spectrogram 

0.12 

1.78 

-7 

Table  VI.1:  Mean  absolute  errors  for  frequency  slope  and  offset.  100  trials,  SNR=10dB, 
performance  breaking  points  for  all  time-frequency  methods  considered. 
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3.b)  Error  Evaluation 

Recall  the  received  signal  was  sampled  at  a  rate  of  512  [samples/pulse  period] 
while  the  image  size  was  set  to  256x256  bins  in  order  to  speed  up  the  computations. 
However,  the  errors  presented  were  corrected  to  correspond  to  an  image  with  size 
512x256,  where  the  time  and  frequency  axes  have  512  and  256  bins  respectively.  Next, 
the  goal  is  to  investigate  how  the  estimated  errors  relate  to  the  frequency  error  obtained 
for  the  frequency  expression  given  in  equation  (VI-8). 

Recall  that  the  normalized  frequency  expression  is  given  by: 

fn  =/wo+«-«> 

and  the  slope  angle  ■&  measured  in  the  image  is: 

a-N-2Nf 

tan(0)  = - L,  (VI-12) 

where  a  is  the  frequency  slope  parameter,  N  is  the  length  of  the  signal,  and  Nf  and 
Nt  represent  the  number  of  bins  in  the  frequency  and  time  axis  respectively.  In  our  case 
Nf=N/2  and  Nt=N.  Thus,  equation  (VI- 12)  simplifies  to: 

tan($)  =  a-  N  .  (VI- 13) 

Recall  that 

a  f  k 

s[n]  =  cos(2n(fn0 -n  +  —  -n1)),  where  fn0=- j~,  a= — -  ,  n=l,..,512 

2  Js  fs 

Let’s  assume  that  the  linear  chirp  has  frequency  slope  k0  and  that  %  is  the 
corresponding  angle,  which  leads  to: 
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(VI- 14) 


k0  _  _  tan(r?0 ) 

—  ~ao  ~ 


N 


Now,  let’s  assume  we  estimate  an  angle  equal  to  tiest  instead  of  t?0due  to  errors  in 


the  estimation  process.  Thus, 


$est  =  t?0+At?, 


(VI-15) 


where  Ai?  is  the  angle  error.  As  a  result,  the  resulting  estimated  frequency  slope 
kesl  expression  is  given  by: 

2  tan(t?0  +  At?) 


*«»=/,  *««=/, 


N 


_  J  s 


fs  tan(t?0 )  +  tan(  At?) 


N  1  -  tan(t?0 )  •  tan(At?) 

ko-f,2-N  +  f,*-  tan(At?) 
fs2  ■  N -k0  ■  N2  ■  tan(At?) 


(VI- 16) 


Thus,  the  final  frequency  slope  error  is: 


Ak=kn  -k^,  = 


(fc02-V2+//)-tan(At?) 


0  est  -  -t2  4._/a.on  r2 


k0-N*  •  tan(At?)  —  //  -V 


(VI-17) 


The  slope  frequency  error  Ak  is  expressed  in  terms  of  the  angle  error,  but  also 
depends  on  the  sampling  frequency,  and  the  length  of  the  signal  N. 

If  we  assume  again  that  f0  is  the  initial  offset  value  corresponding  to  equation  (VI- 
8),  then  the  normalized  offset  frequency  is  derived  from  equation  (VI-1 1)  asfn0=fo/fs.  The 
valued  measured  is  expressed  in  bins  and  is  related  to  fno  by: 

fb0=fn0-2-Nf  (VI-18) 

Assume  the  estimated  initial  offset  value  fbest  contains  an  error  Afb.  Thus, 


54 


f best  f  bO  A/i b  » 


(VI- 19) 


and  the  estimated  value  for  /o  according  to  equations  (VI-11)  and  (VI-8)  is  given 
by: 

/<w  =  /,  •  fb0+Afb  =f0+fs- ■  (VI-20) 

0est  s  2Nf  2  Nf 

When  Nf=N/2,  the  resulting  error  is: 

A/o  =  /o-/o„,=-/,'^L-  (VI-21) 

Equations  (VI- 17)  and  (VI-21)  provide  the  relations  between  errors  measured  in  the  time- 
frequency  image  representation  of  size  N/2xN,  where  N  is  the  length  of  the  signal,  and 
the  corresponding  errors  in  the  initial  time-varying  frequency  expression  given  in 
Equation  (VI-8).  Note  that  Afb  is  measured  in  bins. 

3.c)  Multi-pulse  processing  results 

Up  to  this  point,  all  results  were  derived  by  extracting  the  frequency  modulation 
parameters  from  one  noisy  radar  pulse.  Performances  improve  when  using  more  than  one 
pulse.  Assume  we  isolate  five  radar  pulses,  where  each  pulse  contains  the  same 
transmitted  signal  at  a  given  SNR  level.  Note  that  the  signal  information  gets  mapped  to 
the  same  location  of  the  time-frequency  plane,  while  the  noise  contributions  scatter  to 
different  location  in  each  trial.  Thus,  averaging  multiple  time-frequency  images  improves 
the  quality  of  the  signal  information. 
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Thus,  averaging  the  contribution  of  five  pulses  leads  to  the  results  shown  in 
Figure  VI.  12  obtained  for  the  Smoothed  Pseudo  Wigner-Ville  transformation.  One 
hundred  realizations  per  SNR  level  are  used,  and  SNR  levels  varied  from  -20  to  OdB. 
Note  that  there  is  no  need  to  consider  higher  SNR  levels  as  the  error  curves  already  flatten 
for  SNR=0dB.  Results  show  a  significant  improvement  over  using  one  pulse  only. 
Further  improvements  may  be  obtained  by  considering  a  larger  number  of  pulses. 
However,  this  will  result  in  a  direct  increase  of  the  computation  time. 


Figure  VI.12:  Slope  and  offset  error  for  the  SPWV  distribution;  five  realizations. 

3.d)  Application  to  the  AN/SPS-40B  radar 

We  apply  the  above  results  to  a  specific  radar  using  pulse  compression  readily 
available  to  us:  the  SPS-40B  radar.  The  Radar  Set  AN/SPS-40B  is  used  for  the  detection 
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ranging  and  tracking  of  air  targets  at  long  range  in  American  and  foreign  navy  services 
[24].  In  this  pulse  mode,  the  pulse  length  is  T=60psec.  The  transmitted  waveform  inside 
the  pulse  is  a  linear  modulated  down-chirp  with  a  bandwidth  of  2MHz  centered  at  a 
frequency  that  can  be  varied  manually  by  the  operator  from  402.5  to  447.5  MHz. 

Let’s  assume  that  the  radar  operates  at  437.5  MHz,  i.e.,  the  chirp  starts  with  a 
frequency  438.5MHz  and  decays  linearly  to  436.5  MHz.  Our  ESM  (Electronic  Support 
Measures)  receives  a  series  of  noisy  pulses  from  this  radar.  Further,  assume  we  can 
isolate  one  pulse  perfectly  and  that  the  estimated  received  signal  is  centered  at  438MHz 
with  a  bandwidth  equal  to  2MHz.  Given  that  the  duration  of  the  pulse  is  60psec  and  if  we 
use  512  [samples/pulse  duration],  the  sampling  frequency  should  be  set  at: 

fs  =  ~7ZMHz  =  8-53 MHz  .  (VI-22) 

60 

This  sampling  frequency  does  not  satisfy  the  Nyquist  rate.  Thus,  we  can  either 
increase  the  sampling  frequency  or  heterodyne  to  a  lower  center  frequency.  The  first 
option  has  the  drawback  of  increasing  the  number  of  samples  to  deal  with,  and  the 
computational  load  of  the  estimation  schemes  considered.  Heterodyning  the  signal  down 
to  the  baseband  can  be  accomplished  by  multiplying  the  received  signal  with  a  cosine 
function,  and  using  a  lowpass  filter  to  extract  the  information.  Assume  the  heterodyning 
frequency  is  selected  as  440MHz.  Thus,  the  resulting  chirp  signal  is  an  up-chirp  at  the 
frequency  440-437.5=2.5MHz  with  bandwidth  equal  to  2MHz  after  heterodyning  and 
filtering.  The  theoretical  time-frequency  representation  is  shown  in  Figure  VLB.  Thus, 
the  instantaneous  frequency  of  the  heterodyned  signal  is  given  by: 
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fit)  =  1.5-106  +  3.33-1010  -  t . 


(VI-23) 


The  heterodyned  signal  may  be  sampled  with  a  frequency  fs— 8.53Mhz, 
corresponding  to  512  samples/pulse  duration. 

Assume  we  use  the  SPWV  distribution  and  the  chirp  parameter  estimation 
procedure  described  above  in  Section  V-B.  Figure  VI.  12  shows  that  the  mean  slope  error 
is  0.09°  and  the  mean  offset  error  is  1  bin  when  the  SNR  is  equal  to  5dB.  Using  equations 
(VI-17)  and  (VI-21)  leads  to: 

IAkl=2.6108  and  IAfol=1.67104 , 

meaning  that  the  normalized  errors  for  the  frequency  slope  and  offset  are  equal  to 

_ =  0  008  or  0.8%  and  =0,0104  or  1.04%  respectively  for  the 

3.33. 1010  1.2-106 

heterodyned  signal. 

Note  that  the  normalized  slope  error  remains  the  same  for  the  original  received 
signal,  as  the  signal  bandwidth  has  not  been  changed  by  the  heterodyning  process,  while 
the  normalized  frequency  offset  error  becomes  equal  to  1.67104/438.5106,  i.e.,  0.0038%. 
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Figure  VI.13:  Instantaneous  frequency  for  the  heterodyned  signal;  Normalized  frequency 

/,=  I- 

4.  Comparisons 

Other  schemes  extracting  the  modulation  parameters  have  been  proposed  in  the 
literature  with  their  own  advantages  and  drawbacks  [25-33].  In  this  section,  we  consider 
the  scheme  proposed  by  Peleg  and  Porat  which  estimates  the  parameters  of  a  complex 
linear  FM  scheme  from  a  finite  number  of  noisy  observations 

x[ri\  =  s[n]  +  w[n],n  =  1,...,  N,  where  w[n]  is  white  noise  and  s[n]  is  a  linear  chirp  signal 
defined  as: 

s[n]  =  cos(2 n(f„0n  +  ^n2)),n  =  1 

which  leads  to  the  following  normalized  time-varying  frequency  for  the  chirp: 
fn  =  fno  +an,n  = 

4.a)  Introduction 
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Peleg  and  Porat  proposed  a  fast  estimation  scheme  based  on  two  FFT  operations 
and  two  one-dimensional  searches  of  the  resulting  FFT  quantities  [25].  The  scheme  uses 
the  fact  that  the  maximum  likelihood  estimates  of  the  parameters  of  a  complex  linear  FM 
signal  can  be  derived  by  the  two-dimensional  maximization  of  the  maximum  likelihood 
function  L(a,fn0  )  defined  as: 


L(a,fn0) 


X  exP(-^—  +  /«0n) 

n= 1  * 


Peleg  and  Porat  proposed  to  replace  the  two-dimensional  maximization  operation 
by  two  successive  one-dimensional  maximization  searches  to  obtain  a  sub-optimal 
solution  obtained  with  the  following  steps: 

1)  Compute  the  discrete  form  of  the  ambiguity  function  DAF(x,  CD,n0)  defined  as: 

N-n0 

DAF(x,(d,n0)  =  5>n+T*>-'a\ 

n-\ 

The  parameter  n0  is  chosen  equal  to  N/2,  as  this  value  minimizes  the  mean  square  error  of 
the  slope  parameter  a  as  a  function  of  nQ! N,  as  shown  by  Peleg  and  Porat  [25]. 


2)  Find  the  frequency  co^  which  maximizes  the  magnitude  of  DAF(x,  Ct),n0), leading  to 
- 


a  =  - 


nn 


3)  Define: 

zn=xne(-fan2\n  =  l,...,N. 

4)  Compute  the  parameter  fn0  which  maximizes  L(co,  a)  with  respect  to  co : 

|2 


Uoi.a)  =  \DFT(z„o>)\  = 


1=1 
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Peleg  and  Porat  showed  that  their  scheme  approaches  CRB  values  for  SNR  values  above 
5dB.  They  also  indicated  that  the  SNR  level  should  satisfy  the  following  constraint  for  a 
given  number  of  data  points: 

SNR  >-^(50  +  V2500  +  50  V ) 

where  N  represents  the  data  length,  for  the  results  to  be  meaningful  [25]. 

4.b)  Simulations 

We  generated  one  hundred  linear  chirps  of  length  512  samples  with  random  start 
and  stop  normalized  frequencies  in  the  range  [0,  .5]  for  a  given  SNR  level.  The  SNR  was 
varied  in  the  -lOdB  to  lOdB  in  steps  of  2dB,  as  done  previously  in  Section  VI.C.3,  and 
the  chirp  parameters  estimated  following  the  scheme  proposed  by  Peleg  and  Porat. 
Figure  VI.  14  presents  the  results  obtained.  We  compared  the  results  obtained  for  our  best 
three  transform  types  (Pseudo  Wigner-Ville,  reassigned  pseudo  Wigner-Ville,  and 
reassigned  smoothed  Pseudo  Wigner-Ville  distributions)  to  those  obtained  by  Porat  and 
Peleg.  Results  show  that  the  Porat  and  Peleg  scheme: 

1)  starts  to  break  down  around  -7dB,  and 

2)  is  not  as  robust  to  noise  degradations  as  those  based  on  combined  TFR/Radon 
transforms  as  the  estimation  errors  are  significantly  larger  for  SNRs  below  - 
5dB. 
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Figure  VI.14:  Slope  and  offset  errors  obtained  with  Peleg  and  Porat’  scheme;  linear  chirp 
signals  [25];  100  trials  per  SNR  level;  left  plot:  slope  error;  right  plot:  offset  error. 
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VII.  HYPERBOLIC  CHIRP  PARAMETERS  ESTIMATION 


This  section  considers  the  estimation  of  hyperbolic  chirp  parameters.  As  before, 
the  starting  point  is  the  time-frequency  image  representation  of  the  information. 
However,  generalizations  of  the  Radon  transform  to  hyperbolic  line  types  were  not  used. 
Instead,  we  consider  an  iterative  procedure  to  estimate  the  chirp  parameters. 

Note  that  we  assume  we  know  the  type  of  modulation  transmitted,  as  we  do  not 
address  issues  related  to  distinguishing  between  linear  and  hyperbolic  modulation  types. 
Such  classification  issues  are  left  for  future  work. 


A.  SIGNAL  GENERATION 
1.  Introduction 

Assume  that  we  can  isolate  individual  pulses  of  duration  t=512  samples.  Thus,  a 
received  signal  with  hyperbolic  modulation  frequency  is  given  by: 

x(t)  =  cos(27T(Aln(f  +  t0)  +  B- 1 )) ,  (VI-1) 

where 

/«=  — +  B.  (VE-2) 

t  +  t0 


The  analytical  signal  s(t)  obtained  from  x(t)  with  a  Hilbert  transform  is  given  by: 


x _  ej-2n(Aln(t+t0)+Bt) 


(vn-3) 


The  corresponding  discrete  signal  \[n]  is  given  by: 


;2?r(/»ln(— +^)+B— )  j2x(A-ln(n+nn)-A-]n(fs)+—n) 

x[n]  =  e  fs  fs  fs  —  e  fs  t 


(vn-4) 
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or  finally 


x[n]  =  K  ■  e^aHn+n .H*«)  f  (Vn-5) 

where: 

a  =  A,  b  =  —  ,  n0  =  fs-t0,  and  K  =  g-/'2®4111^  ,  (VII-6) 

f s 

The  normalized  instantaneous  frequency  expression  is  given  by: 

fn[n]  =  -^-+b.  (VD-7) 

n  +  n0 

Recall  that  we  need  to  select  the  parameters  a,  b  and  no  such  that  the  range  of  the 
normalized  frequency  fn[n]  is  between  0  and  0.5  to  prevent  aliasing,  which  leads  to  the 
following  ranges: 

ae[0,oo],  noe[°’H,  be  [-**>,' 0.5]  (VD-8) 

A  valid  selection  for  a,  b  and  n0  is  obtained  by  selecting  one  parameter  in  its 
allowed  range,  and  then  the  other  two  so  that  they  also  fall  in  their  allowed  ranges.  In  the 
simulation,  we  first  select  a  value  for  b  from  a  uniform  distribution  in  the  region  [0,  0.5]. 
The  starting  frequency  at  n=0,fn(0)  ,is  selected  from  a  uniform  distribution  in  the  region 
[b,  0.5].  The  final  frequency  fn(N),  for  n=N,  where  N  is  the  signal  length,  can  be  selected 
in  the  range  [b,fn(0)].  Defining  b,  fn(0),  and  fn(N)  leads  to  the  following  values  for  a  and 
no'. 

a  (/n(0)-frHfr  +  /n( o))-*  (vn_9) 

(fn(o)-fnm) 
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(vn-io) 


(b  +  fn(N))-N 

"°  (/„(0)-/„(AT))' 

Note  that  this  random  selection  can  sometimes  lead  to  a  frequency  equation  with  a 
very  steep  slope  at  the  beginning,  as  shown  in  Figure  VII.  1.  Such  a  behavior  is 
undesirable  because  none  of  the  time-frequency  representations  gives  good  resolution  in 
the  area  of  the  steep  slope,  and  this  modulation  type  is  not  typical  in  radar  applications. 


Theoretical  Time-Frequency  diagram 


Figure  VII.l:  Hyperbolic  modulation;  normalized  frequency  for  a=10,  n0=25,  b=0.05. 

2.  Hyperbolic  Line  Parameter  Constraints 

We  wish  to  restrict  our  random  selection  of  the  signal  parameters  to  those  leading 
to  chirps  without  steep  slopes,  and  need  to  automate  the  process  so  that  we  can  run 
multiple  trials  to  study  the  scheme  robustness  to  noise  degradations.  Thus,  we  define  a 
figure  of  merit  that  describes  the  amount  of  curvature  present  in  a  given  hyperbolic  line. 
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The  parameter  selected  to  characterize  the  curvature  of  the  chirp  is  the  distance  d  defined 
as  the  maximum  distance  between  any  point  of  the  curve  and  its  cord. 

Note  the  parameter  b  shifts  the  hyperbolic  chirp  up  or  down  without  affecting  the 
shape,  thereby  its  curvature.  Therefore,  we  assume  that  b=0  for  simplicity.  In  such  a 

case,  the  chord  associated  to  the  hyperbolic  chirp  fn[n]  =  — - — is  a  line  passing  through 

n  +  n0 

the  points  (0,—)  and  (N, — - — )  with  equation: 
n0  N  +  Hq 

a  a  a 

_ *o  N  +  n0  n0 

n- 0  N- 0 

The  resulting  chord  equation  is  given  by: 

a-n  +  riQ  •  (N  +  nQ)  ■  f  -  a  ■  (N  +  Uq)  =  0,  (VII- 1 1) 

which  leads  to: 

/= - — - n  +  — .  (Vn-12) 

n0-(N  +  n0)  n0 

Recall  that  the  gradient  of  a  curve  at  a  given  point  is  a  line  that  passes  through  that 
point  and  has  a  slope  equal  to  the  derivative  of  the  curve  at  that  point.  Thus,  the  gradient 
at  any  point  of  the  hyperbolic  line  is: 

A  =  f'=  ~~-y  (VH-13) 

(n  +  n0) 
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The  distance  d  is  maximum  at  the  point  n=< \  where  its  gradient  is  parallel  to  the 

chord.  As  a  result,  the  gradient  slope  at  the  point  n=%  is  equal  to  the  slope  of  the  chord. 

Thus  the  position  of  the  point  n=( \  can  be  estimated  from  (VII- 1 1)  and  (VII- 13)  as: 

—  a  _  -a 
(£  +  n0)2  nQ-(N  +  n0)  ’ 


which  yields  the  coordinates  of  the  point  E,  as: 


(Vn-14) 


m= 


(vn-i5) 


Recall  the  distance  d  between  a  point  (xi,yi)  and  a  line  with  equation 
aj  •  x  +  a2  ■  y  +  a3  =0  is  given  by  the  equation: 
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(vn-i6) 


d  = 


I  aj  •  jtj  +  a2  ■  yj  +  a3  I 


4 


2  2 
<*\  +<*2 


Thus,  the  distance  between  the  chord  and  the  point  (£,f(£»  can  be  computed  using 
equations  (VII- 11),  (VII- 14),  (VII- 15),  and  (VII- 16),  which  leads  to: 

\a-£+n0-(N  +  n0)-f(Z)-a-(N  +  n0)\ 

d  —  i  —  , 

^Ja2  +n$(N  +  n0)2 


and: 


a 


(V  +  2-n0)-2--y/n0(N  +  no) 


ija2 +n$(N +  nQ): 


(vn-i7) 


The  distance  d  can  be  viewed  as  a  figure  of  merit  for  the  curvature  of  the 
hyperbolic  line.  A  high  value  of  d  means  the  curvature  is  high  and  the  time-frequency 
representation  near  the  time  origin  is  poor.  However,  very  small  values  of  d  represent 
cases  for  which  the  hyperbolic  curve  is  almost  indistinguishable  from  a  straight  line. 
Thus,  we  restricted  the  chirp  signals  generated  so  that  the  distance  d  is  in  the  region 
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to  avoid  such  cases. 


Next,  additive  zero-mean  white  Gaussian  noise  is  added  to  generate  the  noisy 
chirp  with  a  specific  SNR  level. 


B.  FEATURE  EXTRACTION 
1.  Introduction 

The  signal  time-frequency  representation  can  be  obtained  with  any  of  the  energy 
distributions  discussed  earlier  in  Section  V.  However,  note  that  the  atomic 
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decompositions  are  not  as  well  suited  as  the  energy  distributions  for  the  task  at  hand,  as 
they  do  not  describe  the  hyperbolic  line  curvature  accurately.  Therefore,  we  only  consider 
energy  distributions  in  this  chapter. 

At  this  point  the  task  is  to  extract  the  three  unknowns  parameters  ( a,b,no )  for  a 
given  time-frequency  representation.  The  basic  Radon  transform  can  no  longer  be 
applied,  as  it  is  defined  for  straight  lines  only.  The  Radon  transform  was  extended  to 
detect  functions  of  arbitrary  shape  [17],  however  the  computational  load  is  significantly 
higher. 
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Figure  VII.3:  Extraction  of  the  instantaneous  normalized  frequency  expression  from  the 
time-frequency  image.  Before  and  after  median  smoothing,  L=5. 


69 


The  method  considered  here  approaches  the  problem  quite  differently.  It  is  an 
adaptive  procedure  which  fixes  one  of  the  three  unknown  parameters  at  each  iteration, 
and  estimates  the  other  two.  The  resulting  scheme  is  presented  next. 

2.  Method  Description 

The  instantaneous  frequency  expression  for  the  hyperbolic  chirp  is  extracted  from 
the  two-dimensional  image  by  selecting  the  peak  values  obtained  at  each  time  bin  of  the 
image.  The  result  is  a  vector  containing  the  frequency  values  for  each  time  bin,  as  shown 
in  Figure  VII.3. 

Note  that  the  2-D  instantaneous  frequency  approximation  is  very  accurate  at  high 
SNR  levels,  and  degrades  as  the  SNR  level  decreases.  In  noisy  environments,  the  pixel 
with  the  highest  energy  obtained  at  a  given  time  bin  may  be  an  outlier,  resulting  in  spikes 
in  the  instantaneous  frequency  estimate,  as  illustrated  in  the  middle  plot  of  Figure  VII.  3. 
Such  outliers  can  be  smoothed  out  with  a  median  filter.  Simulations  showed  that  a 
median  filter  of  length  5  smoothed  out  potential  “spikes”  without  loss  of  resolution. 

We  selected  the  three  energy  transformations  leading  to  the  best  time-frequency 
image  quality  for  the  linear  chirp  case,  and  restricted  our  hyperbolic  chirp  analysis  to 
those.  The  distributions  selected  are: 

-  Pseudo  Wigner-Ville  distribution. 

Reassigned  Pseudo  Wigner-Ville  distribution, 

-  Reassigned  Smoothed  Pseudo  Wigner-Ville  distribution. 

We  set  the  dimension  of  the  image  at  512x512  bins  to  increase  the  image 
resolution  and  reduce  quantization  errors.  In  addition,  simulations  showed  that  these  three 
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energy  distributions  do  not  necessarily  perform  well  at  the  beginning  and  end  of  the 
image.  As  a  result,  we  only  consider  the  image  from  time  bin  60  to  time  bin  450. 

The  second  step  in  the  proposed  scheme  approximates  the  instantaneous 
frequency  adaptively  with  a  hyperbolic  line  by  minimizing  the  squared  error  between  the 
information  contained  in  the  image  and  a  theoretical  hyperbolic  curve  expression.  Note 
that  the  problem  to  be  solved  is  non-linear,  due  to  the  specific  frequency  law  to  be 
approximated,  as  we  estimate  the  parameters  a,  b  and  no  given  in  equation  (VII-7).  We 
first  tried  to  solve  the  problem  using  a  classical  nonlinear  least  square  iteration  scheme 
provided  with  the  function  “lsqnonlin”  from  the  MATLAB  optimization  toolbox  [15]. 
Simulations  showed  that  the  algorithm  converged  to  different  local  minimum,  depending 
on  the  initial  values  selected.  However,  this  problem  can  be  resolved  using  a  two-step 
procedure  as  follows. 

Assume  we  wish  to  approximate  the  hyperbolic  curve  given  in  equation  (VII-7) 
with  a  function  of  the  form: 

y(n)  = — — — .  (VH-18) 

n  +  n0 

If  we  assume  our  estimation  values  obtained  from  the  image  to  be  equal  to  yest(n), 
n=l,. . .N,  then  we  wish  to  find  a  and  no  so  that: 

yest(n) - —  =  0,  for  n=l...N.  (VH-19) 

n  +  n0 

The  above  set  of  equation  forms  a  linear  system  of  N  equation  with  two 
unknowns,  which  can  be  solved  using  a  least-square  method.  Next,  assume  we  have  an 
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estimate  of  the  parameter  b,  best,  which  contribution  is  subtracted  from  the  frequency 
equation  given  in  (VIT-7),  resulting  in: 

(vn-20) 

n  +  Uq 

Thus,  the  problem  becomes  to  fit  the  data  f[n]  by  finding  the  parameters  (a,n0) 
which  best  fit  a  curve  of  the  form  a/(n+n0)  in  a  least  squares  sense.  The  set  of  estimated 
parameters  has  a  mean-square  error  e,.  At  this  point,  the  problem  becomes  to  update  the 
parameter  b,  and  re-estimate  corresponding  values  for  a  and  n0  so  that  the  error  function 
expressed  as  a  function  of  b  is  convex  with  a  strong  minimum.  The  location  of  the 
minimum  value  for  the  error  function  is  obtained  for  best  and  the  best  estimated  values  of 
a  and  n0. 

Even  though  we  could  not  formally  prove  that  the  shape  of  the  error  function  as  a 
function  of  b  is  convex,  we  observed  the  same  type  of  convex  shape  for  the  error  function 
over  150  randomly  generated  hyperbolic  chirps.  Figure  VR4  plots  the  normalized  errors 
and  the  mean  values  for  a,  b,  and  no. 
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Realizations 

Figure  VII.4:  Normalized  errors  obtained  for  a,  b,  n0 ;  150  randomly  generated  hyperbolic 
chirps. 

The  mean  and  standard  deviation  of  the  error  is: 


Mean 

Std  Deviation 

a 

0.0361% 

0.000287 

b 

0.1092% 

0.0013 

n0 

0.0262% 

0.000207 

Note  that  the  mean  error  for  b  is  slightly  higher  than  that  of  the  other  two 
unknowns.  This  is  to  be  expected,  as  b  is  restricted  to  the  range  [0,  0.5].  Therefore,  very 
small  error  values  may  correspond  to  large  normalized  errors. 
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We  adaptively  estimate  the  value  for  b,  by  taking  advantage  of  the  convex  shape 
of  the  error  function.  First,  we  restrict  the  search  to  a  specific  region  of  b.  Next,  we  select 
five  values  of  b  equally  spread  within  that  range,  and  compute  the  other  two  parameters 
and  the  resulting  mean-square  error.  Next,  we  restrict  the  b  range  to  that  including  the 
minimum  location  by  using  the  information  provided  by  the  mean-square  errors,  and 
repeat  the  process.  This  iterative  process  can  be  performed  forever,  as  we  have  no 
knowledge  of  the  minimum  mean  square  error.  In  practice,  we  restricted  the  number  of 
iterations  to  10,  as  the  values  of  b  were  restricted  to  small  range  [0,0.5]  in  our  simulations 
to  keep  the  computational  load  under  control.  Theoretically,  the  range  of  b  can  set  as 
large  as  we  want.  The  algorithm  can  converge  in  any  area  of  b  but  it  will  require  a  larger 
number  of  iterations  to  preserve  the  same  accuracy  as  the  range  of  b  expands,  resulting  in 
a  computational  load  increase. 

C.  SIMULATION  RESULTS 

A  few  comments  are  in  order  before  discussing  the  simulation  results: 

1.  The  scheme  considered  above  is  not  as  robust  to  noise  degradations  as  the 
linear  chirp  scheme  described  in  Sections  V  and  VI.  This  is  to  be  expected  as  a  relatively 
clear  and  undistorted  time-frequency  image  is  required  to  extract  the  normalized 
frequency  information. 

2.  The  iterative  scheme  finds  the  set  of  a,  b  and  n0  which  minimize  the  mean- 
square  error.  However,  relatively  large  error  values  in  the  parameter  estimates  may 
correspond  to  small  error  in  the  actual  hyperbolic  curves.  Figure  VII.5  shows  true  and 
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estimated  hyperbolic  curves.  The  true  parameter  values  are:  a- 65,  6=0.1  and  «o= 240, 
while  the  estimated  parameter  values  are:  <3=84,  6=0.08  and  na= 294.  The  normalized 
errors  are  29%,  20%  and  22%  for  a,  b  and  no  respectively,  even  though  the  two  curves  are 
hardly  different. 


Figure  VII.5:  Two  hyperbolic  lines;  Solid  line:  true  hyperbolic  curve  (a=65,  b=0.1  and 
n0=240,  dotted  line:  estimated  hyperbolic  curve  (a=84,  b=0.08  and  n0=294). 


Hyperbolic  chirps  were  generated  by  randomly  selecting  the  parameters  a,  b,  and 
no  in  the  allowed  ranges  mentioned  earlier.  Next,  additive  white  Gaussian  noise  was 
added  to  generate  SNR  levels  between  0  and  20dB,  in  steps  of  2dB.  One  hundred 
realizations  were  generated  for  a  given  SNR  level.  Figures  VII.6  to  VTI.8  plot  the  mean 
and  the  standard  deviation  for  the  normalized  errors  for  (a,  b,  n0)  as  a  function  of  the  SNR 
level,  where  the  normalized  error  is  defined  as: 

I  true  value  -estimated  value  I 

norm  error  = - — - -100  % 

true  value 
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Note  that  the  above  definition  may  lead  to  large  errors  when  the  parameter  true 
values  are  very  close  to  zero,  due  to  the  denominator.  Thus,  we  restricted  our  simulations 
to  cases  where  b  is  in  the  range  [0.025,  0.5].  The  other  two  parameters  a  and  n0  are 
selected  using  the  method  described  in  section  (VII- A). 


Figure  VII.6:  Hyperbolic  chirp;  normalized  errors  for  the  SPWV  distribution. 

Figures  VH6  to  VII.  8  show  the  normalized  error  mean-square  and  corresponding 
standard  deviation  for  (a,b,n0)  obtained  using  the  Smooth  Pseudo  Wigner-Ville,  the 
Reassigned  Pseudo  Wigner-Ville,  and  the  Reassigned  Smoothed  Pseudo  Wigner-Ville 
distribution  respectively. 
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Results  show  that  the  normalized  errors  decrease  to  zero  as  the  SNR  level 
increases.  They  also  show  that  the  SPWV  is  the  scheme  most  robust  to  noise  degradations 
out  of  the  three  considered. 


SNR  (db) 

Figure  VII.7:  Hyperbolic  chirp;  normalized  errors  for  the  RPWY  distribution. 

Results  also  show  that  the  reassigned  methods  perform  better  than  the  SPWV  for 
high  SNRs.  This  is  to  be  expected,  as  they  provide  a  more  focused  image  by  finding  the 
center  of  gravity  of  the  energy  distribution  for  each  time  instance,  resulting  in  a  better 
image  quality.  However,  the  reassignment  process  worsens  the  image  quality,  as  the  noise 
level  increases,  resulting  in  the  estimation  process  breaking  down.  Simulations  show  that 
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the  RPWV  performs  slightly  better  than  the  RSPWV  does,  especially  for  SNR’s  in  the 
range  of  2  to  5dB. 
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Figure  VII.8:  Hyperbolic  chirp;  normalized  errors  for  the  RSPWV  distribution. 


The  continuous  parameters  A,  B  and  to  in  equation  (VII-2)  are  related  to  the 
parameters  of  the  discrete  signal  expression  via  equation  (V1I-6).  Thus,  the  normalized 
errors  estimated  for  the  discrete  values  a,  b  and  no  are  identical  to  those  obtained  with  the 
continuous  parameters. 
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VIII.  CONCLUSIONS 


This  study  investigated  the  detection  of  pulse  start  and  stop  times  of  noisy  pulsed 
chirp  signals  and  the  estimation  of  chirp  modulation  parameters  for  two  specific 
modulation  schemes;  linear  and  hyperbolic. 

The  robustness  of  three  TCF-based  schemes  and  an  envelope  detection  algorithm 
in  noisy  environments  were  compared.  Results  showed  that  none  of  the  pulse  detection 
schemes  considered  in  this  work  to  be  clearly  better  than  the  others  for  the  range  of  SNR 
levels  considered,  and  the  specific  selection  to  be  a  function  of  the  desirable  characteristic 
(letitbePFA  or  Pd)  to  be  optimized. 

The  estimation  of  the  modulation  parameters  was  approached  from  an  image 
processing  point  of  view.  We  compared  eleven  different  time-frequency/scale 
transformations  and  investigated  their  robustness  to  represent  noisy  chirps.  We  compared 
the  quality  of  the  estimated  modulation  parameters  obtained  for  linear  chirps  when 
applying  a  Radon-based  transformation  to  the  time-frequency/scale  images.  Results  show 
time-frequency  transformations  to  lead  to  better  focused  images  when  dealing  with  noisy 
chirp  signals,  and  to  better  estimation  of  the  modulation  parameters  than  wavelet-based 
decompositions  do.  Specifically,  the  Smoothed  Pseudo  Wigner-Ville  distribution  had  the 
best  performance  in  low  SNR  environments,  while  the  reassignment  scheme  improved 
the  performances  for  high  SNRs.  The  combination  of  the  energy  distributions  and  the 
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Radon  transform  lead  to  small  detection  errors  and  almost  constant  performances  for 
SNRs  as  low  as  -7dB. 

Results  also  showed  that  the  estimation  performances  increase  when  using 
multiple  pulses.  Simulations  showed  that  the  parameter  detection  SNR  threshold  dropped 
to  -lOdB  when  processing  five  pulses  simultaneously.  However,  the  computation  time 
increases  significantly. 

Hyperbolic  chirp  parameters  were  extracted  from  the  time-frequency  image  using 
an  iterative  procedure.  However,  this  scheme  requires  a  very  good  estimate  of  the 
instantaneous  frequency  expression  as  an  initial  estimate,  thereby  restricting  its 
application  to  medium  to  high  SNR  levels. 

This  study  is  restricted  to  two  specific  types  of  modulation;  linear  and  hyperbolic. 
Additional  work  would  be  needed  to  extend  this  type  of  approach  to  additional 
modulation  types.  In  addition,  we  assumed  the  modulation  type  to  be  known  a-priori. 
However,  modulation  identification  issues  need  to  be  investigated  to  allow  for  an 
automated  scheme  to  be  implemented.  A  significant  amount  of  work  dealing  with 
modulation  identification  schemes  is  available  in  the  literature,  each  with  their  own 
advantages  and  drawbacks  [45-47].  Further  work  is  needed  to  investigate  how  a 
modulation  identification  scheme  can  be  integrated  with  the  work  discussed  here. 
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